rm(list=ls())

load("TSCS_data_draws.RData")

no <- paste("X",1:1000,sep = "")

# OLS ---------------------------------------------------------------------

linear_list <- list()

for (variable in no) {
  
  formula = as.formula(paste0("downturn ~ lag(downturn,1:2) +
      lag(antipluralism_gov_seat_share,1)*",
                              lag(variable,1),
                              "+ lag(polarization,1) + 
      lag(clientelism_gov_seat_share,1) + 
      lag(region_citizen_support,1) + 
      lag(region_antipluralism_gov_seat_share,1) + 
      lag(region_liberal_democracy,1) + 
      lag(democratic_stock,1) + 
      lag(GDP_capita,1) + 
      lag(life_expectancy,1)"))
  
  out <- plm(formula,
             data = panel_draws,
             model = "within",
             effect = "twoways")
  

  linear_list[[variable]] <- out
  
}

### plot results

coeffs_lm <- vector()
se_lm <- vector()

for (i in no) {
  coeffs_lm[i] <- linear_list[[paste0(i)]][["coefficients"]][[paste0("lag(antipluralism_gov_seat_share, 1):",i)]]
}

vcm_lm <- lapply(linear_list, function(x) plm::vcovBK(x, cluster=c("group")))
rse_lm <- lapply(vcm_lm, function(x) sqrt(diag(x)))

for (i in no) {
  se_lm[i] <- rse_lm[[paste0(i)]][[paste0("lag(antipluralism_gov_seat_share, 1):",i)]]
}

coeff_df <- as.data.frame(cbind(coeffs_lm,se_lm))
coeff_df %>% mutate(margin_of_error = se_lm*1.65,
                    CI_lower = coeffs_lm - margin_of_error,
                    CI_upper = coeffs_lm + margin_of_error) %>%
  arrange(coeffs_lm) %>%
  mutate(model = row_number())->
  coeff_df

coeff_df %>% 
  mutate(sig = ifelse(CI_upper < 0,1,0)) %>%
  group_by(sig) %>%
  summarise(n = n())

# 844 significant, 156 insignificant

pdf(file = "fig_e1a.pdf",width = 5,height = 9)

ggplot(coeff_df,aes(x=coeffs_lm,y=model))  + 
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), colour="grey") +
  theme() + xlab("Coefficient") + ylab("Model")  +
  geom_vline(xintercept = 0) + geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

dev.off()

# GMM ---------------------------------------------------------------------

gmm_list <- list()

for (variable in no) {

  panel_draws$interaction <- lag(panel_draws$antipluralism_gov_seat_share, 1) * lag(panel_draws[, variable], 1)
  
  formula_str <- paste("downturn ~ lag(downturn, 1:2) + lag(antipluralism_gov_seat_share, 1) + lag(citizen_support, 1) +",
                       "interaction + lag(polarization, 1) + lag(clientelism_gov_seat_share, 1) +",
                       "lag(region_citizen_support, 1) + lag(region_antipluralism_gov_seat_share, 1) +",
                       "lag(region_liberal_democracy, 1) + lag(democratic_stock, 1) + lag(GDP_capita, 1) +",
                       "lag(life_expectancy, 1) | lag(downturn, 3:5)", collapse = " ")
  
  formula <- as.formula(formula_str)
  
  out <- pgmm(formula,
              data = panel_draws,
              effect = "individual",
              model = "onestep",
              transformation = "ld",
              indexes = c("country", "year"),
              robust = TRUE)
  
  gmm_list[[variable]] <- out
}

### plot results

coeffs_gmm <- vector()
se_gmm <- vector()

for (i in no) {
  coeffs_gmm[i] <- gmm_list[[paste0(i)]][["coefficients"]][["interaction"]]
}

vcm_gmm <- lapply(gmm_list, function(x) plm::vcovHC(x, cluster=c("group")))
rse_gmm <- lapply(vcm_gmm, function(x) sqrt(diag(x)))

for (i in no) {
  se_gmm[i] <- rse_gmm[[paste0(i)]][["interaction"]]
}

coeff_df <- as.data.frame(cbind(coeffs_gmm,se_gmm))
coeff_df %>% mutate(margin_of_error = se_gmm*1.65,
                    CI_lower = coeffs_gmm - margin_of_error,
                    CI_upper = coeffs_gmm + margin_of_error) %>%
  arrange(coeffs_gmm) %>%
  mutate(model = row_number())->
  coeff_df

pdf(file = "fig_e1b.pdf",width = 5,height = 9)

ggplot(coeff_df,aes(x=coeffs_gmm,y=model))  + 
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), colour="grey") +
  theme() + xlab("Coefficient") + ylab("Model")  +
  geom_vline(xintercept = 0) + geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

dev.off()

